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SUMMARY 

Incompatible meshes, i.e. meshes that physically must have a common boundary, 
but do not necessarily have coincident grid points, can arise in the course of a 
finite element analysis. For example, two substructures may have been developed at 
different times for different purposes and it Dec ones necessary to interconnect the 
two models. A technique that uses only multipoint constraints, i.e. MPC cards (or 
MPCS cards in substructuring), is presented. Since the method uses only MPC's, the 
procedure may apply at any stage in an analysis; no prior plaining or special data 
is necessary. 


INTRODUCTION 

When two separate finite element solid meshes are combined with a common bound- 
ary, the degrees of freedom on the common surface of one of the models must be 
eliminated. The elimination process must preserve the compatibility of the joined 
surface when the model is subjected tc external loads which cause the interface sur- 
face to deform. The use of MPC's can provide this capability to first order, i.e. 
linear in terms of grid point displacements. 

The MPC equations presented are based on a lineal' surface spline. A three 
dimensional surface spline is a mathematical tool used to find a function, f{u,v,w) 
for all points (x,y,z) when *( u(x,y,z), v(x,y,z), w(x,y,z) is known for a discrete 
set of points, (x^y^zj) a..J (u,v,w) are displacements in the (x,y,z) directions. 

The three dimensional surface spline given here is based on the shape function of 
the linear displacement tetrahedra. This shape function will produce compatible 
surface deformations for models using linear elements. 

The example presented implements this method by assuming that one surface is 
composed of a s^t of imaginary triangles formed by grid points on one of the mating 
surfaces. The vertices of a triangle form a unique surface described by (x^,y^,z^), 
where i = 1,2,3. Grid points of the surface to oe interconnected that lie within 
the boundary of the triangle are MPC f ed to the displacement of the vertices, 
( u i* v i» w i)» based on the spline functions f( u(x^,yj[,Zi ), v(xi,yi,Zi), v(xi,yi,Zi) ). 
Results of a practical example in substructure analysis are presented. A computer 
program that automatically writes the MPC cards is included in the Appendix. 

DEFINITION OF THE PROBLEM 

Concider a planar triangle in space defined by three points; A, B and C, as 
shown in Figure 1. 

The vector [A} represents the distance from the origin to point A; similarly 
for points B, C and 0. The subject of this paper is the determination of the 
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( 1 ) 


relationship between displacements of points A, B, C and 0 when 

[dA > = [A' - / } = (u a , v*, w 4 ) 

[dB}= [B' - B } = (u b -. " b , v b) 

[dC } = (C' - C } = (uc, y«, v c ) 

[d0>= [0* - 0 } = (u°, v°, w°) 

such that the point 0 # lies in the plane described by the triangle in space defined 
by the three points A’, B f and C*. 

In the application described, the points A, B and C lie on the surface of a 
region described try a group of finite elements; the point ”0" is a grid (or node) in 
the same or in a different group of finite elements defining a region that has a 
physically congruent undeforraed boundary that the points A, E, C and 0 lie on. 

In the notation used here, { ] is a row vector and { ) is a scalar quantity. 

A matrix multiplication, {][][} results in a scalar, as indicated by the left 
and right elements; similarly, [ > { ] would indicate a matrix result. In the text 

a scalar will often be written without the braces. 

SELECTION OF THE SPLINE FUNCTION 

Some rule must be selected that relates the displacement of any point in the 
plane of the triangle with the motion of the three reference points A, B and C, i.e. 
a spline function. The shape function of a linear displacement tetrahedra is used 
for this purpose. However, in order to use this shape function, a tetrahedra must 
be constructed from the reference triangle. This is easily accomplished by forming 
a cross product of the vectors [V^Jand [V^}; 

[yapj = [yah} x [V“} (2) 

where 

[Vij} = [VJ} - [vi} 

The result of equation ( 2 ) is zhe formation of a reference point, P, which forms the 
fourth V2rtex of the reference tetrahedra. Clearly, if [v a °} and [V* 0 } are unique, 
then point P will be unique. 

The shape function for a linear displacement tetrahedra is [Reference 1 , Section 

5.12] 

‘1 x a y* z a l F D x 

1 x b y b z b D 2 ( 3 ) 

= 1 x c y° z c D 3 

.1 xP yP zPJ LDj* 

or 

[u 1 } = [T][D> 

where 

u* = the displacement of the i tb vertex in the x coordinate direction 

Two additional sets of equations, similar to equation ( 3 ), relate v and w, the 
y and z displacement, to the undetermined coefficients, [D}. The displacement of a 
point M 0 lf is, in the x direction 

{ u°} = Di + D 2 X° + D3Y 0 + D^Z 0 (k) 
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or from the solution of equation (3) 

tU | — \-L X J ^ Jl^j iuj 

where 

[s] = [T]-l 

Similar equations exist for v and w. 


It is useful to write equation (U) in the following form for the displacements 
in the three coordinate directions: 

u°(x°, y°, z°) = /: u a + B u^ + C u c + P u P 

v°(x°, y°, z°)=Av a + Bv t> + Cv c + PvP (5) 

w°(x°, y°, z°) = Aw a + Bw* ) + Cw c + PwP 

where 


A = S n + S 21 x° + S 31 yO + S kl ?° 
B - S ^2 + S 22 x° + S 32 y° + Si| 2 z° 
C = S 13 ♦ S 2 3 x ° + $33y° + Si, 32 ° 
P - S ^4 + S 2 i|X° + ^i^y 0 + SI 4 I 4 Z 0 

and 


u* = displacement of point (x*, y*, z*) 
v* = displacement of point (x*, y\ z ] ) 
w* = displacement of point (x*, y*, z i ) 


in X direction 
in Y direction 
in Z direction 


Now equation (5) is in the form of a multipoint constraint (MPC) relationship 
[Reference 1, Section 3.5] • The coefficients A, B, C and P are constants that depend 
only on the constrained geometry. Equation (5) cannot be used directly, however, 
since point P does not actually exist. Point P must be removal from the set of 
equations. 

The motion of point P, in terms of the motion cf point A and the rotation vec- 
tor, [a, 6, y}, of the triangle is [Reference 1, Section 3-5]: 



K) 

r° 

z a P 

-y ap "| fai 


Up = 

I V s 1 + 

- z ap 

0 

xap B 

(6) 

LwP) 

|v*l 

L y* p 

-x a P 

0 J Ly) 


[uP> = 

[u a } + 

[H] [a } 





where 

a, 6, and y = rotation about x, y, z axes. 


CALCULATION OF THE ROTATION VECTOR 


It is necessary to determine the relative angular displacement of the un- 
strained and strained triangles in terms of the displacements of points A, B and C. 
Equation (6) can be written in matrix form for the three points A, B and C and 
the rotation vector (a, 8, y) as: 
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0 z ab — y^ b 1 


Q 1 

u ab \ 

& 

0 

0) 

to 

1 


6 = 

vab | 

yab -xab o 



yab | 

0 z^ -yhc 


ubc / 

C 


v bc l 

ybc - x bc 0 


w bc f 

o 2 ca -yea 


u ca 

-z ca 0 x ca 


vca 

yea - x ca o 




or, symbolically 
[R] [a } = [ Au) 

vher 

[ uar l fu r - u a ) 

v&r = I v r - vM 

var) [_w r - v*) 

Row a unique solution to equation (7) exists so long as the triangles (A, B, C) 
and (A 1 , B f , C') are unique. Physically this must be the case, otherwise a zero 
strain would result. The solution is [Reference 2]: 

U> = ([S]T (R])~l [RJT [ Au > ( 8 ) 

Equation (8) gives the rotation vector, for small angles, when the displacements 
are known at the triangle vertices. 

Equation (6), when [a> is determined from equation (6) becomes: 

[uP} = [u a } + [Q][Au> (9) 

where 

[Q] = [h](rTr)-i [r]T 

Equation (9)* when combined with equation (5) provides the MPC relationship 
that is a function of the displacements of the original four points of the tetra- 
hedra. A, B, C and P. 


Equation (9) written in detail is: 



The displacement of point P can be determined from the displacements u a , u° 
and u c in accordance with equation (9)* i.e. 
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( 10 ) 


uP = (1 + Qjj - Qn)u a + (Qn - Qilj)u^ + (Qii* - Qiy)u c 

+ t'Qlg - Qx2)v® + (Q12 “ 0 x 5 )v tl + (Q15 - QiqJv 0 

+ (Qly - Qi^w 4 + (Q13 - Qi6)w b + (Qig - Q^Jw 0 

and similarly for vP and vP. 

SPLINE EQUATION FINAL FORM 

Substituting equation (10) into the MPC equation ( 5 ) gives the final form of 
the spline equations; 

u(x°, y°, z°) = 

[A + P(l+Q 17 - Q n )]u a + [J^P(Q n - Qm) ] u b + [O+PCQm - Q 17 }] u<= 

+ P(Qig ” ^12^ v 3, + p (Ql2 “ Qlg) v b + P(Ql5 ~ Qig) v° 

+ P(Qi 9 - Q13) w® + P(Q 13 - Q l6 ) w b + P(Q l6 - Q 19 ) wC 

v(x°, y°, z°) = 

P(02 7 - Q2l) u a + P(Q2 x - Q2 J*) u b + P(Q 2 U - Q27) u c 

+[A + P(l+Q2g _ ^ 2 ^ 1 ya+ [B+P(Q2 2 “ ^ 25^1 ' /b + [C+P(^5 - Qgg)] v° 

P(029 ~ 023^ + P(^23 - 026^ v b + P(Q26 ~ Q29) w 0 

w(x°, y°, z°) = 

P( 0 37 “ 0 3 i) u a + P( 0 3 i - Q3I4 ) u b + P( Q3I4. - 0 37 ) u c 

+ P( 0 3 g “ 0 3 2 ^ v 3 + P(Q 3 2 “ Q35) v b + ^(§35 ” Q 3 g) v° 

+[A +P(l+Q39 - Q 33 ) Iw 3 - + [B+P(Q33 - Q36)] w b + [C+P(Q 3 g - Q39)] w° 


A SIMPLE NUMERICAL EXAMPLE 


Consider a triangle with vertices at A=(l,0,0), B=(0,1,0) and C=(0,0,l). The 
point, 0, to be removed is at 0=(l/3, 1/3, 1/3). This point lies in the plane of 
the triangle ABC with the distance to the origin of l/yH . Values of the various 
scalars, vectors, and matrices used in calculating the MPC coefficients are given 
below. 


[yap } = 
[T] = 


( 1 . 0 , 

1.0, 

1.0) 


'1.0 

1.0 

0.0 

0.0' 

1.0 

0.0 

1.0 

0.0 

1.0 

0.0 

0.0 

1.0 

. 1.0 

2.0 

1.0 

1.0 
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•v o w > 


[ S ] = 

' .6667 

.3333 

.3333 

-.3333 


.3333 

-.3333 

-.3333 

.3333 


-.6667 

.6667 

-.3333 

.3333 


L— .6667 

-.3333 

.6667 

.3333 

A = 

i/3 




B = 

1/3 




C = 

1/3 





0 



([ R t ][ R]) -1 = f .2778 -.0556 -. 0556 ' 
-.0556 .2778 -.0556 

-.0556 -.0556 .2778 


([ rJ t [ r])-i [ R] t = 


.0556 

.0556 

.222 

-.m - 

.2778 - 

.2778 

.0556 

.2222 

. 0556 ' 

.0556 

.0556 

.222 

.222 

.0556 

.0556 - 

.2778 -.1111 - 

.2778 

-.2778 

-.2778 

-.111 

.222 

.0556 

.0556 

-0556 

, 2222 

. 0556 . 

' .3333 

.3333 

• 3333 

.0000 

.0000 

.0000 

-.3333 -.3333 

-. 3333 ' 

-.3333 

-.3333 

-.3333 

.3333 

.3333 

.3333 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.3333 

-.3333 

-.3333 

.3333 

.3333 

.3333 
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and finally , 

u(l/3, 1/3, 1/3) = 

• 3333u a + .3333u b + .3333u c 


+ .0 

v* + 

.0 v b 

+ 

.0 v° 

+ .0 

w® + 

.0 v b 

+ 

.0 wC 

v( 1/3, 1/3, 1/3) 

= 



.0 

U a + 

.0 u b 

+ 

.0 u c 

+ .3333V® + 

. 3333v b 

+ 

• 3333v c 

+ .0 

v a + 

.0 v b 

+ 

.0 v° 

v(l/3, 1/3, 1/3) 

= 



.0 

u a + 

.0 u b 

+ 

.0 u c 

+ .0 

va + 

.0 v b 

+ 

.0 v c 

+ . 3333w® + 

. 3333w b 

+ 

. 3333w° 


If the displacements of points A, B, and C each move along the x, y and z axes, 
respectively, as shown in Figure 2, then 

u° = .3333, vQ = .3333, w° = .3333 

and point 0 f = (.6667, .6667, .6667). It is useful to no^e that since P = 0 the 
MPC equations (ll) reduce to: 

u(x°, y°, z°) = A u a + B u* 3 + C u c 

v(x°, y°, z°) = A v 3 - + B v^ + C v c 

Thir. situation apparently results when the point to be removed, 0, lies in the 
plane of the triangle formed by points A, B and C. The equations derived are such 
that the strain perpendicular to the triangle ABC are null. Points lying in the 
plane of the triangle have non-null strains only in the plane of the triangle. 

A PRACTICAL APPLICATION 

The technique presented in this paper was developed during the course of a 
finite element analysis of a bolted, flange-type coupling used in marine risers 
[3,10. Figure 3 shows the coupling components in the bolt-up condition. The 
coupling is modeled as three separate substructures; the Pin, Box, and Bolt, as 
shown in Figures U, 5, and 6. Advantage is taken of the cyclic symmetry of the 
coupling, hence a 22.5° pie section is modeled with the appropriate boundary 
conditions. The finite element mesh of the Pin in the area of the countersink used 
very small elements in order to properly calculate stress concentrations. The 
corresponding area of the Bolt, at the Pin-Bolt interface, need not be modeled with 
a mesh this fine. 
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A considerable cavings in computer run time is achieved if the interface of the 
Pin and Bolt can be made consistent with the boundary conditions but allow a lower 
density mesh on the Bolt. Grid points on the Pin that were not coincident with 
grid points on the Bolt were removed. The results of a stress analysis of the Pin 
are shown in Figure 7. These results show that, in the region of the grids that 
were removed by MPCS cards, a smooth and consistent stress pattern results. 
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FIGURE 1 - SCHEMATIC REPRESENTATION OF AN UNSTRAINED TRIANGLE <A,B,C) 
AND A STRAINED TRIANGLE (A',B',C') 



FIGURE 2 - EXAMPLE PROBLEM 



FI6URE 3 - HHF RISER CONNECTOR 





FIGURE 4 - HIDDEN LINE VIEW OF HflF PUL 22.5* Pit SECTION 
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FIGURE 5 - HIDDEN LINE VIEW OF HNF BOX, 22.5* PIE SECTION 
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FIGURE 6 - HIDDEN LINE VIEW OF HHF BOLL 22.5* PIE SECTION 
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FIGURE 7 - VON MISES EQUIVALENT STRESS 
CONTOURS MAPPED ONTO A HIDDEN L T NE PLOT 
OF HPT PIN 
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APPENDIX 


This Appendix gives the source code of a FORTRAN computer 
program that reads a control file (BDYMPC. INP) , the NASTRAN 
grid card file ( BDYMPC. BDF) , and produces MPC cards for inclusion 
to the NASTRAN Bulk Data File. Figure A-l shows the control 
file for the sample problem. The first line gives the number 
of points, NG, to be MPC'ed and the three GRID ID's. The next 
NG lines contain the GRID ID's. The group may be repeated as many 
times as necessary. Figure A-2 shows the GRID cards and Figure A-3 
shows a listing of MPC cards generated for the sample problem. 

The source listing is given in Figure A-4. 
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1 21 22 23 

25 


FIGURE A-l SAMPLE PROBLEM COMTBOL FILE 
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GRIE 

GRID 

GRID 

GRID 


21 

22 

23 

25 


0 1 .# 0.0 0.0 

0 0.0 1.0 0.0 

0 0.0 0.0 1.0 

0 0.33333 0.33333 0.33333 


FIGURE A-2 SAMPLE PROBLEM GRID CARD INPUT 
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$ DEBUG 

SNOFLOATCALLS 

PROGRAM BO TMPC 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OMIT 

ii 

* . INP 

OMIT 

12 

*.BOF 

OMIT 

13 

* .MPC 

OMIT 

14 

*.DIA 


CONTROL INPUT DATA, CORNER GRIO ID'S 
NO. AND ID'S OF INTERNAL GRIDS 

NASTRAM GRID CARDS OF FEM 

OUTPUT OF NASTRAM MPC CARDS 

DIAGNOSTIC AND CHECK OUTPUT 


IMPLICIT REAL *8 (A-H.O-Z) 

DIMENSION AP(3) .0(3,9) ,R (9,3) ,T (4,4 ) ,S (4,4) ,RT (3,9) ,RTR(3,3) 
fc,RTRI(3,3),H(3,3),RTRIRT(3,9) 

CHARACTER DOM* 8 

COMMON /ZMFCR/XMPC ( 18 ) , IDG ( 3 ) , LIME , X ( 3 ) 

OPrN ( II , FILE* * BOTMPC. INP ', STATUS" ' OLD * ) 

OPEN ( 12, FILE" 1 BOTMPC. BDF ' , STATUS" 'OLD * ) 

OPEN (13, FILE" 'BOTMPC. MPC* , STATUS ■' NEW • ) 

OPEN ( 14 , FI LE" * BDTM PC . DI A ' , STATUS " 1 NEW • ) 

NCR"0 

LINE"1 

XMPC ( 1 ) "-1 . ■ 

DO 108 1-1,4 
100 T(I,1)"1.0 

C 

C START LOOP FOR EACH REGION 
C 


105 READ (11,5,END"1000)NG,IDG 
5 FORMAT (BN, 41 8) 

WRITE ( 14, 1)NG, IDG 

1 FORMAT (IX,* NO. OF INTERNAL GRIDS • • ,18 
A, ' CORNER GRID ID"S ',318) 

NCF-0 


C FIND CORNER GRID ID'S 

C 

107 REWIND 12 

110 READ (12 , 2,END"2000 ) DOM, ID , IOUM,X 
2 FORMAT (BM,A8,2I 8, 3E8.0) 

WRITE (14, 4) DOM,ID,IOOM,X 
4 FORMAT (IX, AS, 218, 3E12. 5) 

IF (ID.NE.IOG(NCF-M)) GO TO 110 
WRITE (14, '(IX)') 

NCF-NCF+1 
DO 120 I "1,3 
120 T(HCF,I+1)-X(X) 

IF(NCF.EQ.3) GOTO 130 
GOTO 107 
130 CONTINUE 


C 


FIGURE A-4 COMPUTER PROGRAM SOURCE CODE LISTING 
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MFC 


312 

2S 

1 

•l.MI 

21 

1 

• 333 

♦M 

1 

♦M 

1 


22 

1 

.333 

23 

1 

• 333 

♦M 

2 

♦M 

2 


21 

2 

.IN 

22 

2 

.888 

♦H 

3 

♦M 

3 


23 

2 

.999 

21 

3 

.888 

♦M 

4 

♦M 

4 


22 

3 

.999 

23 

3 

• 888 



MFC 


312 

25 

2 

• 1.999 

21 

1 

.088 

♦M 

S 

♦M 

5 


22 

1 

.999 

23 

1 

.088 

♦M 

6 

♦M 

6 


21 

2 

.333 

22 

2 

.333 

♦M 

7 

♦M 

7 


23 

2 

.333 

21 

3 

.888 

♦M 

8 

♦M 

8 


22 

3 

.Iff 

23 

3 

.088 



MFC 


312 

25 

3 

-l.fff 

21 

1 

.888 

♦N 

9 

♦M 

9 


22 

1 

.HI 

23 

1 

.888 

♦M 

18 

♦M 

18 


21 

2 

.HI 

22 

2 

.088 

♦M 

11 

♦M 

11 


23 

2 

.Iff 

21 

3 

.333 

♦M 

12 

♦M 

12 


22 

3 

.333 

23 

3 

.333 




FIGURE A 3 SAMPLE PROBLEM MPC CARD OUTPUT 
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FORM CROSS PRODUCT A8 X a5 FOR VECTOR AP 

XAB»T(2,2)-T(1,2) 

XAC«T(3,2)-T(1,2) 

XBCT(3,2)-T(2,2) 

YABT(2,3)-T(1,3> 

TAC-T(3,3)-TU»3) 

TBC«T(3.3)-T<2,3) 

*AB-T<2,4)-T(l,4> 

ZAC-T{3,4)-T(1»4) 

ZBC«T(3,4l-T(2»4) 

AP ( 1 ) «<YAB*ZAC-ZAB*YAC 
AP(2)»ZAB*XAC-XAB*ZAC 
AP ( 3 ) »XAB* Y AC-t AB* XAC 

POINT P ADDED TO T MATRIX 

T(4,2)-T (1,2)+AP<1) 

T(4 r 3)*T(l,3)+AP(2) 

T(4,4)»T (l r 4)+AP(3) 

CALL MINV4 (T,S) 

FORM R AND H MATRICES 
NOTE THAT ZCA*-ZAC, ETC. 

R{1,1)« 9.0 
:.(I,2)*+ZAB 
R<1,3)*-YAB 
R (2#I)»-ZAB 
R(2,2)» 9.0 
R(2 r 3)>+XAB 
R (3» 1)»+YA8 
R(3,2)*-XAB 
R(3,3)» 9.0 
R(4 f 1)* 0.0 
R(4,2)«*ZBC 
R<4,3)— YBC 
R(5,l)— Z8C 
R(5#2)« 0.0 
R(5*3)*+XBC 
R(6, 1)«+TBC 
R(6,2)» XBC 
R<6,3)« 9.0 
R(7,l)« 0.0 
R(7,2)— ZAC 
R(7,3)«+YAC 
R (8, i)«+ZAC 
R(8,2)« 0.0 
R(8,3)— XAC 
R(9,1)—YAC 
R(9,2)«+XAC 
R<9,3)» 0.0 
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H(l,l>* M 
H(2«2)* 0.0 
H(3»3)* 0.0 
H(l,2)* AP(3) 

H(l,3)— AP(2) 

H(2,3)* AP(1> 

H(2,l)*-H(l,2) 

H(3,l)**H(l,3) 

H(3,2)*-H(2,3) 

00 135 11*1,3 
00 135 31*1,9 
RT (II, Jl) -R (J1,U) 

135 CONTINUE 

T 

FORM R * R - RTR 

DO 136 11*1,3 
DO 136 Jl*l,3 
RTR(I1,J1)*0.0 
DO 136 Kl*l,9 

RTR(I1,J1)-RTR(I1,J1)+RT(I1,K1)*R(K1,J1) 

136 CONTINUE 

CALL MINV3 (R?R,RTRI ) 

FORM RTRI*RT 

DO 137 11*1,3 
DO 137 Jl*l,9 
RTRIRT (II, Jl J*0.0 
DO 137 Kl-1,3 

137 RTRIRT (II, J1)*RTRIRT (II, J1)+RTRI (I1,R1)*RT(K1,J1) 
C 

DO 138 11-1,3 
DO 138 Jl*l,9 
Q(I1,J1)*0.0 
DO 138 f *1,3 

138 Q (11, Jl)*Q(Il, J1)+H(I1,K1)*RTRIRT (Kl, Jl) 

C 

WRITE (14,3) ( (R (IL, JL) ,JL-1,3) ,IL*1,9) 

..,({ T(IL,JL),JL*1,4),XL*1,4) 

4,(( RTR (IL,JL) , JL*1,3) ,IL*1,3) 

4,(( RTRI (1L,JL) ,JL*1,3) ,IL*1,3) 
k ,AP 

t, ( (RTRIRT (IL,JL) ,JL*1,9) ,IL*1,3) 

4,(( Q (IL , JL) , JL*1 , 9) , IL*1 , 3 ) 

3 FORMAT (IX,’ MATRIX R ’,/9(/3F8.4) 

6 ,//lX, ' MATRIX T ’ ,/4(/4F8.4) 

6, //IX,’ MATRIX RTR ',/3(/3F8.4) 
k ,//lX, * MATRIX RTRI * ,/3(/3F8.4) 

6,//lX,' VECTOR AP ' ,/3F8. 4 
k ,//lX, ' MATRIX RTRIRT ' ,/3(/9F8.4) 
k ,//lX, ' MATRIX Q ' ,/3 (/9F8.4) ) 
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START LOOP FOR BACH GRID TO BE MPC'ED - GRID ID IS THE KEY 

DO 200 IMG-1,MG 
READ (11,5, EHD a 3000 ) IDR 
NCR«NCR+1 
REMIND 12 

READ (12 , 2 , END-2000 ) DON , ID , IDUM , X 
ir (ID. ME. IDR) GOTO 140 

CALC OF A,B,C,P 

ALL ARE FUNCTIONS OF GRID COORDINATE TO BE REMOVED, X (1) ,X (2) ,X (3) 

A*S (1,1) 

B»S (1,2) 

C»S (1,3) 

P-S(l,4) 

DO 150 I>1,3 
A«A+S (X+1,1)*X(I) 

B»B+S (I+1,2)*X(1) 

C*C+S (I+1,3)*X(I) 

P*P+S (X+1,4)*X(I) 

150 CONTINUE 

WRITE (14,6) A,B,C,P 

6 FORMAT (IX, /IX,’ A« ' ,F10.4,' B- ' ,F10.4, ' C- ' ,F10.4 
6,' P- ' ,F10.4,/1X ) 

CALC COEF FOR DIRECTION 1 AND WRTIE MFC CARDS ON UNIT 13 

XMPC(2)« P*(l+Q(l,7)-Q(l,l) )+A 
XMPC (3 )* PM 0(1,1)-Q(1,4) )+B 
XMPC(4 )■ P*( Q(l,4)-Q(l,7) )*C 

XMPC (5)* PM Q(l,8)-Q(l,2) ) 

XMPC (6)* PM 0(1,2)-Q(1,5) ) 

XMPC(7)» PM Q(l,5)-0(i,8) ) 

XMPC (8)* PM Q(l,9)-Q(l,3) ) 

XMPC(9)« PM Q(1,3)~Q (1,6) ) 

XMPC(10)-PM Q (1,6)-Q (1,9) ) 

CALL MPCR (1, IDR) 

LINE>LXMB+4 

CALC COEF FOR DIRECTION 2 AND WRITE MPC CAROS ON UNIT 13 

XMPC ( 2 ) » PM 0 (2,7)-Q (2,1) ) 

XMPC (3)* PM 0 (2,1}-C (2,4) ) 

XMPC(4)« P*( Q'2,4)-Q(2,7) ) 

XMPC (5)" PM 1+0 (2, 8) -Q(2, 2) )+A 
XMPC (6)* PM Q (2,2)*Q (2,5) )+B 
XMPC(7)» PM Q (2,5)-Q (2,8) )+C 
XMPC (8) • PM Q (2,9)-Q (2,3) ) 

XMPC (9) ■ PM Q(2,3)-0(2,6) ) 

XMPC (10 )"P* ( Q(2,6)-Q (2,9) ) 

CALL MPCR (2, IDR) 

L?NE a LINE+4 
C 
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CALC COBP FOR DIR1CTI0M 3 AMD WRITE MFC CARDS OM UNIT 13 

XMFC(2)> P*( Q(3,7)»Q(3,1) ) 

X MFC(3)« PM Q(3,l)-Q(3,4) ) 

XMFC(4)- PM Q(3,4)»Q(3 7) ) 

XMFC(S)- PM Q(3,8)-Q<3,2) ) 

XMPC(C)- PM '}(3,2)‘*Q(3,S) ) 

XMFC (7)* PM Q(3,5)«Q<3,8) ) 

XMFC (8) * PM1+Q(3,9)-Q(3,3) )+A 
XMFC(9)- PM Q(3,3)-Q<3,«) )+Bl 
XHPC(10)-PM Q(3,6)-Q(3,9) ) +C 
CALL MPCR (3# IDR) 

LINE-LINE+4 

200 COMTIKOE 
GOTO 105 

1000 STOP 

2000 WRITB (14 * * ) ’ ERROR NO. 2000* 

STOP 

3000 WRITE (14,*) ’ ERROR NO. 3000' 

STOP 

END 


SUBROUTINE MINV4 (C,A) 

A-INV(C) C IS NOT DESTROYED. ONE«A*C FOR CHECK 

IMPLICIT REALM (A-H,0-Z) 

DIMENSION ONE (4,4) 

DIMENSION A(4,«),C(4,4),B(4,4) , IPVOT (4) , INDEX (4, 2) r PIVOT (4) 

N-4 

M>0 

OO 999 1*1, N 
DO 999 J*1,N 
999 A(I,J)»C(I,J) 

WRITE (14,1) ( (A(I , J) ,J*1,4) ,1*1,4) 

1 FORMAT (IX, 'INPUT MATRIX TO MINV4 ' ,/lX,4 (/1X,4F12.4) ,/lX ) 

57 DET-1.0D0 

DO 17 J*1,M 
17 IPVOT (J) “0.000 
DO 135 1*1, M 

FOLLOWING 12 STMTS SEARCH FOR PIVOT ELEMENT 

T-0.0D0 
DO 9 J*1,N 

IF (IPVOT (J)-l) 13,9,13 
13 DO 23 K-1,N 

IF ( IPVOT (K)-l) 43,23,81 
43 IF (ABS (T) -ABS (A(J,K) ) ) 83,23,23 

83 IROW*J 
I COL 9 K 
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T-A(J,X) 

23 CONTINUE 
9 CONTI NUB 

IPVOT (ICOL) a IPVOT (ICOL) +1 

FOLLOWING 15 STMTS TO PUT PIVOT BLBMBNT ON THE DIAGONAL 

IP (I ROW* ICOL) 73,199,73 
73 DBT—DBT 

DO 12 L a l,N 
T a A(IROW,L) 

A ( I BOW , L) “A ( ICOL , L) 

12 A(XCOL,L) a T 

IF (M) 199,199,33 
33 DO 2 L a l,N 
T a B(IROW,L) 

B < IROW, L)-B( ICOL, L) 

2 B (ICOL,L) *T 

199 INDEX (I ,l)“IROW 
INDEX (I ,2) a ICOL 
PIVOT (I ) -A (ICOL , ICOL) 

DET«DET*PIVOT(X) 

NEXT 6 STMTS TO DIVIDE PIVOT ROW BY PIVOT ELEMENT 

A(ICOL,ICOL) a 1.9D9 
DO 295 L-1,N 

A(ICOL,L) a A(XCOL,L)/PIVOT(I) 

US CONTINUE 

IF (M) 347,347,66 
66 DO 52 L*1,M 

52 B ( ICOL, L) a B (ICOL, L) /PIVOT (I ) 

NEXT 19 STMTS REDUCE NON PIVOT ROWS 

347 DO 135 LX>:,N 

IF (LI -ICOL) 21,135,21 
21 T a A (LX , ICOL) 

A (LI , ICOL) a 9 . 9D9 
DO 89 L"1,N 

89 A (LX ,L) a A (LI ,L) -A (ICOL,L) *T 

IF(M) 135,135,18 

18 DO 68 L«1,M 

68 B(LI,L)-B(LI,L)-B(ICOL,L)*T 
135 CONTINUE 

NEXT 11 INTERCHANGE COLUMNS 

222 DO 1 I-1,N 
L-N-I+l 

IF (INDEX (L,l) -INDEX (L,2) ) 19,3,19 

19 JROW a INDEX(L,l) 

JCOL a INDEX (L , 2 ) 

DO 549 K*1,N 
T a A (X, JROW) 
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A (K, JRON) -A (K, JCOL) 

A(K,JCOL)*T 
549 CORTXNOB 
3 CONTINUE 

si continue 

WRITE (14,992) C (A(L,M) ,H-1,N) ,L*1,N) ,OBT 
992 FORMAT (/IX, * INVERSE OF MATRIX IH MI MV 4 * ,/lX,4 (/1X,4F8.4) 
*,/lX,' DBT T » MPII2.4) 

C 

C CHECK INVERSE AMO MBITS TO (WIT 14 

C 

DO 128 1*1, M 
DO 128 J*1,M 
OMB (I ,J) *8. 808 
00 128 K-1,N 

OMB *,X , J) -acre (1 , J) +A II , X) *C(K, J) 

128 COMTIMOB 

WRITE <14, 991) M,0NS 

991 FORMAT ( IX ,//lX, 1 CHECK I MV ORDER OF MATRIX « *,I4,//1X 
*,4(/lX,4F8.4)) 

RBTQRM 

BHD 

C 

C 

c 

SUBROOTIHB NINV3(T,S) 

C 

C S*IMV(T) T IS MOT OBSTROTEO. 0MS*S*T FOR CHECK 

C 

IMPLICIT REAL *8 (A-H,0~Z) 

DIHEMSIOM T (3,3) ,S (3,3) ,0MB (3,3) ,B(3) ,C (3) ,A (4,4) 

DATA 0ME/9*8.8/ »B/3*8.8/,C/3*B.B/ ,A/1C*8. 8/ 

C 

C MAKE SORB FIRST PIVOT ELBHBMT IS MOT ZERO 

C 

N*4 

NN-W-1 

WRITE (14, 3) ((T(L,H) ,H*1,NM) ,L*1,MH) 

3 FORMAT (IX,//, IX,' MATRIX T * ,/lX,3 (/1X,3F8. 4) ) 
A(1,1)«1.8D8 
00 18 1*1, MM 
DO 18 J*1,MM 
18 A(I+1, J+1)*T(I,J) 

A(l,l)*l .8D8/A(1,1* 

DO 118 M*1,MH 
K*M*1 

58 00 88 1*1 ,M 

B ( I ) •■ . 808 
DO 88 J*1,H 

88 B(I)*B(I)+A(I,J)*A(J,K) 

0*8.808 
00 78 1*1 ,M 
78 D*0+A(K,I)*B(I) 

D*-D+A(K,K) 

A(K,K) *1.8/0 


270 




DO 84 I-1,M 

* * A(I,Kj— B(I)*A(K,K) 

DO M J*1 »M 
C(J)-4.4D4 
00 94 I-l r M 

94 C(J)*C(J)+A(K,I)»A(I. J) 

00 144 J-1,M 
144 A(K.J)*-C(J)*A(K,K) 

00 114 I— 1,M 
00 114 J-l.M 

114 A(I,J)-A(I,J)»B(I)*A(K,'l) 

DO 111 1*1, MM 
CO 111 J*1,NM 
111 S(I,J)-A(I+1,J*1) 

WRITE (14,2) ((A(L,M) ,H*1,M) ,L*1,N) ,0 
2 FORMAT (/IX,* INVERSE OP AUGMENTED T * ,/lX,4 (/1X,4PE.4) 

A, /IX, * DET T • ' . 1FE12.4) 

C 

C CHECK INVERSE AMO UNITE TO OMIT 14 
C 

00 124 1*1, NK 
00 124 J*1,NN 
ONE (I ,J) *4.808 
00 124 K*1,NM 

OME(I,J)-ONE(I,J)+S(I,K)*T(K,J) 

124 CONTINUE 

WRITE (14,1) MM, ONE 

1 FORMAT (IX, //IX, ' CHECK INV ORDER OF MATRIX - *,I4,//1X 
4,3(/lX,3P8.4)) 

RETURN 

END 

C 

c 

c 

SUBROUTINE MPCR (ICOR, IIOR) 

IMPLICIT REAL *3 (A-H,0-Z) 

CHARACTER MFC* 8, SPACE *3, PM* 2 

C O M M ON /XMFCR/ XMEC(IS) ,I0G(3) ,LINE,X(3) 

DATA MPC, SPACE, FN/*MFC ' ,'+M'/ 

ISIO-312 

11*1 

12*2 

13-3 

WRITE (14,1) IIOR, ICOR, X, XMFC 
1 FORMAT (IX,* MFC EQUATIONS: GRIO R E MOVE D *, 


Alt,' DIRECTION ' 

,I4,4X,/1X,4P12.4, 


1/1X,F9.4, * OA 

* ,F9.4, ' OB 

* *P9.4, ' 

UC ', 

2/lX,F9.4, * VA 

* ,F9.4, ' VB 

* ,P9.4, ' 

VC ', 

3/lX,P9.4, * NA 

' ,-9.4, ' MB 

' ,F9.4, ' 

WC ') 

Ll-LIME 




L2-LINE+1 




L3-LIME+2 




L4-LXME+3 




WRITE (13,2) MFC, 

ISIO, IIOR, 

ICOR, XMFC(l) 


1 , IOG ( 1 ) , I 1 , XMFC ( 2 ) ,SPACS , PM, Ll , PM, LI , SPACE , IOG ( 2 ) , I l , XMFC ( 3 ) 
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2,IDG(3),U,XMFC(4),SPACE,PM,L2,PM f L2,SPACE,IDG(l) .I2,XMFC(5) 
3 , tOC ( 2 ) , 1 2 , XM PC ( 6 > , S P ACE , PH , L3 , PH , L3 , S P ACE , I DG ( 3 ) , 1 2 , XM PC ( 7 ) 
4 , IDG ( 1 ) , 13 , XMFC (8 ) , SPACE , PM, L4 , PH, L4 , SPACE , IOC (2 ) , 13 , XMPC (9) 
5* IDG (3) *13, XM PC (1C) 

2 FORMAT (1X,A8,3I8,PC.3 

1.218, F8.3,A8,A2,I€,/1X,A?,I6,A8,2I8,F8.3 

2.218, F8.3,A8,A2,I6,/1X,A2,I6,A8,2I8,F8.3 

3.218, F8.3,A8,A2,I6,/1X,A2,I6,A8,2I8,F8.3 

4.218, F8.3,A8,A2,I6,/1X,A2,I6,A8,2I*8,F8.3 

5.218, r8.3) 

BHD 
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